
function [fd1, po, theored1, red90, red95, red99,red90global,red95global, red99global] ...
    = bpl(data, nw, nzeropad, fmax, plot_x_period, logfreq, linlogY, unit)

% Bending Power Law (BPL)
%
% INPUT
%   data:       evenly spaced data, 2 columns, 1st = depth or time value in ascending order; 2nd = signal value
%   nw:         time-bandwidth product of discrete prolate spheroidal sequences, Typical choices for NW are 2, 5/2, 3, 7/2, 4.
%   nzeropad:   zeropadding number, default = 1024. If nzeropad value is smaller
%               than data length, no padding is applied.
%   fmax:       max frequency, for plotting; default = nyquist frequency.
%   plot_x_period: 1 or 0; default: 1 = plot x axis in period domain; 0 = plot x axis in frequency domain.
%   logfreq:    1 or 0; default: 1 = frequency axis in log scale; 0 = frequency axis in linear scale
%   linlogY:    1 or 0; default: 1 = power axis in log scale; 0 = power axis in linear scale
%   unit:       String; default = no unit
%
% OUTPUT
%   fd1:        frequency
%   po:         power
%   theored1:   bending power law
%   red90:      local red noise 90% confidence level
%   red95:      local red noise 95% confidence level
%   red99:      local red noise 99% confidence level
%   red90global:    global red noise 90% confidence level
%   red95global:    global red noise 95% confidence level
%   red99global:    global red noise 99% confidence level
%   
% EXAMPLES
% data(:,1) = (1:1024)';
% data(:,2) = randn(1024,1);
% nw = 2;
% nzeropad = 1000;
% fmax = 1;
% plot_x_period = 0;
% logfreq = 1;
% linlogY = 1;
% unit = 'm';
% [fd1, po, theored1, red90, red95, red99,red90global,red95global, red99global] ...
%    = bpl(data, nw, nzeropad, fmax, plot_x_period, logfreq, linlogY, unit)
%
% Reference:
%   Vaughan, S., Bailey, R., Smith, D., 2011. 
%   Detecting cycles in stratigraphic data: Spectral analysis 
%   in the presence of red noise. Paleoceanography 26.
%
% From Acycle software (acycle.org)
%
%   by Mingsong Li (Peking University)
%   date: June 6, 2022
%
if nargin < 8; unit = ''; end
if nargin < 7; linlogY = 1; end
if nargin < 6; logfreq = 1; end
if nargin < 5; plot_x_period = 1; end

timex = data(:,1);
datax = data(:,2);
dt = median(diff(timex));

if nargin < 4; fmax = 1/( 2 * dt ); end
if nargin < 3; nzeropad = 1024; end
if nargin < 2; nw = 2; end


df = 1/(timex(end)-timex(1));
bw=2*nw*df;


padtimes = nzeropad/length(datax);

if padtimes > 1
    if nw == 1
        [po,w]=pmtm(datax,nw,nzeropad,'DropLastTaper',false);
    else
        [po,w]=pmtm(datax,nw,nzeropad);
    end
else
    if nw == 1
        [po,w]=pmtm(datax,nw,'DropLastTaper',false);
    else
        [po,w]=pmtm(datax,nw);
    end
end
fd1=w/(2*pi*dt);



% power law
N = length(datax);
% With no padding or smoothing applied
% Number of independent Fourier frequencies
Nf = N/2; 

pol = log(po);
pol = real(pol);
fun = @(v,f)(v(1)*f.^(-1*v(2)))./(1 + (f/v(3)).^(-1*v(2)) + (f/v(4)).^(v(5) - v(2)));
v0 = [20,2.1,5e-2,8e-3,-2];%Might need updating, this is a 1st guess to values you will fit
x = lsqcurvefit(fun,v0,fd1(7:end),pol(7:end));%Careful where to start from
theored1 = exp(fun(x,fd1)); % bending power law fit
K = 2*nw -1;
nw2 = 2*(K);

% Chi-square inversed distribution
% local c.l.
red90 = theored1 * chi2inv(0.90,nw2)/nw2;
red95 = theored1 * chi2inv(0.95,nw2)/nw2;
red99 = theored1 * chi2inv(0.99,nw2)/nw2;

% periodogram method global
alpha90 = 0.10/Nf;
alpha95 = 0.05/Nf;
alpha99 = 0.01/Nf;
red90global = theored1 * chi2inv((1-alpha90),nw2)/nw2;
red95global = theored1 * chi2inv((1-alpha95),nw2)/nw2;
red99global = theored1 * chi2inv((1-alpha99),nw2)/nw2;

% figbpl = figure;
% set(gcf,'Color', 'white')
% if plot_x_period
%     pt1 = 1./fd1;
%     plot(pt1(2:end),po(2:end),'k-','LineWidth',1);
%     if strcmp(unit, '')
%         xlabel('Period')
%     else
%         xlabel(['Period (',num2str(unit),')']) 
%     end
%     xlim([1/fmax, pt1(3)]);
%     set(gca, 'XDir','reverse')
%     hold on
%     plot(pt1,red99global,'r-.','LineWidth',1)
%     plot(pt1,red95global,'r--','LineWidth',2)
%     plot(pt1,red90global,'r-','LineWidth',1)
%     plot(pt1,red99,'b-.','LineWidth',1)
%     plot(pt1,red95,'b--','LineWidth',2)
%     plot(pt1,red90,'b-','LineWidth',1)
%     plot(pt1,theored1,'k-','LineWidth',2)
% else
%     plot(fd1(2:end),po(2:end),'k-','LineWidth',1);
%     xlim([0 fmax]);
%     if strcmp(unit, '')
%         xlabel('Frequency')
%     else
%         xlabel(['Frequency (cycles/',num2str(unit),')'])
%     end
%     hold on
%     plot(fd1,red99global,'r-.','LineWidth',1)
%     plot(fd1,red95global,'r--','LineWidth',2)
%     plot(fd1,red90global,'r-','LineWidth',1)
%     plot(fd1,red99,'b-.','LineWidth',1)
%     plot(fd1,red95,'b--','LineWidth',2)
%     plot(fd1,red90,'b-','LineWidth',1)
%     plot(fd1,theored1,'k-','LineWidth',2)
% end
% legend('Power','99% global','95% global','90% global','99% local','95% local','90% local','BPL')
% ylabel('Power')
% title([num2str(nw),'\pi MTM & bending power law','; bw = ',num2str(bw)])
% set(gcf,'Name',[num2str(nw),'\pi MTM & bending power law'])
% set(gca,'XMinorTick','on','YMinorTick','on')
% 
% if linlogY == 1
%     set(gca, 'YScale', 'log')
% else
%     set(gca, 'YScale', 'linear')
% end
% if logfreq == 1
%     set(gca,'xscale','log')
% end